## Copyright (C) 2012 Estela Moura
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## ACO_TSP_final
## Ant Colony Optimization to TSP

## Author: Estela Moura <estela@kubuntu-Inspiron>
## Created: 2012-06-13

function [MelhorCusto, MelhorPercurso] = RascunhoACO_TSP(MatrizCidades, Iteracoes)

% Parâmetros padrões sugeridos pelo autor do programa original
Formigas = size(MatrizCidades,1);
ExpAlpha = 1;
ExpBeta = 9;
Coef_Evap = 0.9;

## Criar matriz distância
  for(i=1:size(MatrizCidades, 1))
    for(j=1:size(MatrizCidades, 1))
      if(i==j)
	MatDistancia(i,j)=0;
      else
	MatDistancia(i,j)=sqrt((MatrizCidades(j,2)-MatrizCidades(i,2))^2 + MatrizCidades(j,3)-MatrizCidades(i,3)^2);
      endif
    endfor
  endfor
  
  d=MatDistancia;
  n=max(size(d));

  m=Formigas;
  t_max=Iteracoes;

  [L_nn, P_nn] = NearestNeighborTSP(d);

  L_best = inf;	
  T_best = 0;

% INITIALIZATION ===========================================================

% pheromone trails
  c = 1 / (n * L_nn);
  tau = ones(n,n) * c;

% place m ants in n nodes
  ant_tours = zeros(m, n+1);
# # # ant_tours(:,1) = randint(m,1,[1,25]); -> esta função "randint" não existe no Octave
# # # 	aqui, cria-se uma matriz mx1 com valores entre 1 e 25, uma solução equivalente em Octave está abaixo:
  ant_tours(:,1)=randi(25, m, 1);

  t = 1;
  while (t <= t_max)

    % CREATE TOURS =============================================================

    for s = 2 : n
      for k = 1 : m
	current_node = ant_tours(k,s-1);
	visited = ant_tours(k,:);
	to_visit = setdiff([1:n],visited);
	c_tv = length(to_visit);
	p = zeros(1,c_tv);
	for i = 1 : c_tv
	  p(i) = (tau(current_node,to_visit(i)))^ExpAlpha * (1/d(current_node,to_visit(i)))^ExpBeta;
	end
	sum_p = sum(p);
	p = p / sum_p;
	for i = 2 : c_tv
	  p(i) = p(i) + p(i-1);
	end	
	r = rand;
	select = to_visit(c_tv);
	for i = 1 : c_tv
	  if (r <= p(i))
	    select = to_visit(i);
	    break;
	  end
	end
	city_to_visit = select;
	ant_tours(k,s) = city_to_visit;
	tau(current_node,city_to_visit) = (1 - Coef_Evap) * tau(current_node,city_to_visit) + c;
      end
    end

% ENCONTRA OS MELHORES =================================================

  ant_tours(:,n+1) = ant_tours(:,1); %fecha o caminho
  L_T = zeros(1,m);
  best_ant = 1;
  for k = 1 : m
    P = ant_tours(k,:);
    L = 0;
    for i = 1 : n
      L = L + d(P(i),P(i+1));
    end
    L_T(k) = L;
    if (L_T(k) < L_T(best_ant))
      best_ant = k;
    end
  end
  L_min = min(L_T);
  T_min = ant_tours(best_ant,:);

% update pheromone trails;	
  for i = 1 : n
    tau(T_min(i),T_min(i+1)) = (1 - Coef_Evap) * tau(T_min(i),T_min(i+1)) + Coef_Evap / L_min;
  end

% COMPLETE =================================================================

  t = t + 1;
  current_cities = ant_tours(:,n);
  ant_tours = zeros(m, n+1);
  ant_tours(:,1) = current_cities;
  if (L_min < L_best)
    L_best = L_min;
    T_best = T_min;
  end

  end % ends while
  
  MelhorCusto=L_best;
  MelhorPercurso=T_best;

endfunction
